Mon. Not. R. Astron. Soc. 000, [Tp7| (2008) Printed 24 September 2009 (MN LMgX style file v2.2) 



A torque formula for non- isothermal Type I planetary 
migration - I. Unsaturated horseshoe drag 

S.-J. Paardekooper 1 *, C. Baruteau 2 , A. Crida 1 ' 3 and W. Kley 3 

1 DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, United Kingdom 
2 UCO/Lick Observatory, UC Santa Cruz, 1156 High Street, Santa Cruz, CA 95064, USA 

3 Institut fur Astronomie & Astrophysik, Universitdt Tubingen, Auf der Morgenstelle 10, 72076 Tubingen, Germany 



Draft version 24 September 2009 



ABSTRACT 

We study the torque on low-mass planets embedded in protoplanetary discs in the 
two-dimensional approximation, incorporating non-isothermal effects. We couple lin- 
ear estimates of the Lindblad (or wave) torque to a simple, but non-linear, model of 
adiabatic corotation torques (or horseshoe drag), resulting in a simple formula that 
governs Type I migration in non-isothermal discs. This formula should apply in op- 
tically thick regions of the disc, where viscous and thermal diffusion act to keep the 
horseshoe drag unsaturated. We check this formula against numerical hydrodynamical 
simulations, using three independent numerical methods, and find good agreement. 

Key words: planetary systems: formation - planets and satellites: formation. 



1 INTRODUCTION 

Planets are thought to form in circumstellar discs around 
young stars. In the core accretion model, gas giant planets 
emerge in the disc through gas accretion onto a previously 
accumulated solid core of a few times the mass of the Earth 
(Mj, Pollack et al. 19961. An alternative scenario involves 
direct fragmentation of the disc ( Boss|1997 1 , which is proba- 



bly only possible in the outermost regions of the disc (Bolcy 



2009) 



In general, objects embedded in protoplanetary discs 
will exchange angular momentum with the disc, which leads 
to a change in their orbital parameters. The nature of this 
interaction depends on the masses of the object and the disc. 
Small bodies, up to a few km in size, on a Keplerian orbit 
will experience a head wind from the gas, since the gas is 
partially supported by pressure and will thus orbit at sub- 



Keplerian velocity ( Weidenschilling 1977). This head wind 



will lead to orbital decay, the time scale of which can be as 
short as a few 100 yrs (Weidenschilling 1977). 

The most massive objects, approximately the mass of 
Jupiter, can tidally truncate the disc, forming a deep annu- 
lar gap around their orbits ( Lin fc Papaloizou||1986a| ). The 
planet, being repelled by both gap edges, is locked inside 
the gap and will slowly accrete onto the central star with 
the rest of the disc (Lin & Papaloizou 1986b I. The mini- 



mum mass for this Type II migration to occur depends on 
the scale height and viscosity of the disc ( Crida et al.| 2006 ) . 
Planets that are not massive enough to open up a gap, 
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but do significantly perturb the disc, can be subject to a 
very rapid mode of migration called Type III when embed- 
ded in a very massive disc ( jMasset fc Papaloizou|2003| ) . The 
mechanism of Type III migration relies on a distortion of 
streamlines in the coorbital region due to a radial flow of 
gas with respect to the planet jPepliriski 2008]) . This radial 
flow of gas can be due to the migration of the planet it- 
self, resulting in a positive feedback with the possibility of a 
runaway process ( jMasset fc Papaloizou|[2003[ ), with migra- 
tion time scales of the order of a few tens of dynamical time 
scales. Sustaining this rapid mode of migration has proved 
to be very difficult ( |Peplinski et al.|[20"08a|b[ ). 

In this paper, we will be concerned with planets that 
do not significantly perturb the disc, which is typically valid 
for objects up to a few Mj . This regime of Type I migration 
was long thought to be the simplest case, since it could be 
treated using a linear analysis. It was shown in |Goldreich fc| 
|Trema inc ( 1979) that the torque exerted on the planet by the 
disc can be decomposed in a wave torque, arising at Lind- 
blad resonances, and a corotation torque, generated at coro- 
tation resonances. This analysis was subsequently refined 



(Artymowicz 1993 Ward 19971, eventually resulting in a 



semi-analytical torque formula for isothermal discs (Tanaka 



|et al.|2002| ). This formula has been confirmed by fully non- 



linear, isothermal, hydrodynamical calculations (Bate et al. 
2003||D'Angelo et al.|2003 l. 



The time scale for Type I migration is inversely propor- 
tional to the mass of the planet and the disc, but is typically 



10" 



yr for a 1 Mq planet embedded in a Minimum Mass 



Solar Nebula ( |Ward|1997| |Tanaka et al.|2002| >. This is wor- 
rying, since the lifetime of the disc is of the order of 10 6 
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yr, making the survival of low-mass planets highly unlikely. 
Planetary synthesis models have great difficulties reproduc- 
ing the observed semi-major axis distribution when includ- 
ing Type I migration, and need to reduce the Type I torque 
from an order of magnitude ( Ida fc Lin|2008 1 to as much as 
a factor of 1000 jAlibert et al.|2005| |Mordasini et al.|2009 |. 

Above results were obtained without considering mag- 
netic effects. It has been shown that including magnetic 
fields, either regular ( |Terquern||2003[ ) , or turbulent ( |Nelson| 



equation of motion take the form 



& Papaloizou 20041, may slow down or even stop Type I 



migration. It is not clear, however, if protoplanetary discs 
are sufficiently ionised throughout to couple effectively to 
the magnetic field. 

There is another important ingredient missing in stan- 
dard models of Type I migration, which is to release the 
isothermal assumption and account for the energy balance in 
a more realistic way. A growing body of studies is dedicated 



to this problem, dealing with high-mass planets (D'Angelo 
et al.||2003| |Klahr fc Kley||2006| >, shadowing-effe cts (|Jang- 



Condell & Sasselov 20051, and opacity jumps (Menou & 



Goodman|2004l). InlPaardekooper & Mellema|(|2006al), it was 



shown through three-dimensional, radiation-hydrodynamic 
simulations that for deeply embedded low-mass planets, 
Type I migration could be qualitatively different from the 
isothermal case. Planets could suddenly move outward as 
well as inward, depending on the local opacity. This result 
was confirmed using two-dimensional simulations with a self- 
consistent heating and cooling balance ( |Kley fe Crida| 2008). 
It was subsequently shown that this effect was due to the 
effect of a radial entropy gradient in the disc on the coro- 



tation torque ( jBaruteau fc Ma ssct 2008a; Paardekooper & 
Mellema||2008[), a nd non-linear in nature ( |Paardekooper fc 



Papaloizou|2008 1 



In this series of papers, we aim at catching the essential 
physics of the non-linear, non-isothermal corotation torque 
in a simple model that can be used to predict the Type I 
migration rate, as a function of radial density and temper- 
ature gradients. In this paper, we consider the unsaturated, 
adiabatic horseshoe drag, combined with a linear estimate 
for the wave torque. Effects of viscous and thermal diffusion 
will be considered in a forthcoming work. We start in section 
[2] with reviewing the basic equations and disc models, and 
describe our numerical methods in section|3] We give a more 
detailed overview of isothermal Type I migration in section 
[4] In section [5] we present a simple model for the torque 
on a low-mass planet in the presence of both entropy and 
vortensity gradients, and subsequently compare this model 
to numerical simulations in section [6] A short discussion is 
given in section [7] and we present our conclusions in section 

El 



2 BASIC EQUATIONS 
2.1 Governing equations 

The basic equations are those of the conservation of mass, 
momentum and energy for a two dimensional disc in a frame 
rotating with angular velocity Q p . We adopt a cylindrical 
polar coordinate system (r, ip) with the origin (r = 0) lo- 
cated at the central mass. The continuity equation and the 



9E 
dt 
and 
Dv 

7h 



= -V ■ (Ev) 



+ 2f2 p k x v 



-Vp - V$ 



(1) 



(2) 



respectively, while, in the adiabatic case, entropy is con- 
served along streamlines: 



Dt 



0, 



(3) 



where 7 is the adiabatic exponent. Above, E denotes the 
surface density, v the velocity, p is the vertically integrated 
pressure, $ is the gravitational potential and k is the unit 
vector in the vertical direction. The convective derivative is 
defined by 



D d „ 

— = h v ■ V. 

Dt dt 



(4) 



In the remainder of this paper, we will refer to s e Jj/E 7 
as the entropy of the fluid. An ideal gas equation of state 
was used, p — R g ET/^i, where R g is the gas constant, n is 
the mean molecular weight and T is the temperature. We 
neglect effects of self-gravity, viscosity and thermal diffusion. 
The potential $ contains terms due to the central mass M* , 



and direct and an indirect term due to the planet (see Nelson 
et al |2000| >. 



2.2 Equilibrium models 



We construct axisymmetric equilibrium models that have 
power law profiles in surface density and temperature, with 
indices —a and — /3 respectively. This means that the initial 
entropy profile is a power law as well, with index — £, where 



£ = /3-(7-l)a. 



(5) 



The angular velocity is Keplerian, with a slight correction 
for the radial pressure gradient to maintain pressure equilib- 
rium. The temperature at the location of the planet is cho- 
sen so that the pressure scale height at the location of the 
planet is H p = hr p , with h <C 1. Typically we use h = 0.05. 
In the absence of self-gravity, the density at the location of 
the planet E p can be chosen arbitrarily. 



2.3 Planet 

The potential of the planet, located at r — r p and tp = ip p , 
is taken to be a softened point mass: 



GM P 



r p — 2rr p cos((^ — ip p ) + b 2 r p 



(6) 



with b the softening parameter. In order to approximately 
account for 3D-effects, b should be comparable to h. Typi- 
cally we use b = OAh. When calculating the torque on the 
planet, we include all disc material. We have checked that 
excluding a fraction of the Hill sphere in the torque calcu- 
lation does not affect the results. Below, we will use q to 
denote the mass ratio M p /M». 
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3 NUMERICAL METHODS 

Equations [l] [2] and [3] are solved on a cylindrical grid, ex- 
tending from r/r p = 0.4 to r/r p = 1.6, and the full 2ir in 
azimuth. The typical resolution amounts to Ar/r p — 0.0013 
and A<p = 0.0025. Due to the small radial extent of the 
horseshoe region, a large radial resolution is required. We 
have checked that taking square cells around the planet's 
location (by doubling the resolution in ip) does not influence 
the results. 

We have used three independent numerical codes: 
RODEO (ROe solver for Dis c Embedded Objects, 
Paardekooper & Mellema 2006b), based on an approxi- 



numerically using outgoing wave boundary conditions (Ko 



mate Riemann solver, FARGO (Fast Advection in Rotating 
Gaseous Objects, [Massct 2000a b), and RH2D (Radiation 
Hydrodynamics in 2 Dimensions, Klcy 1989, 1999). The lat- 
ter two methods are based on the van Leer upwind algo- 
rithm. 

RODEO is based on the general relativistic Roe solver 



outlined in Eulderink & Mellema (19951. It uses station- 



ary extrapolation to integrate gravitational and geometrical 
source terms, and can handle arbitrary coordinate frames. 
Since it is based on a Riemann solver, RODEO is specifically 
designed to handle sharp discontinuities, usually in the con- 
text of shocks. We will see in section [6] that although shocks 
do not play a role for low-mass planets, discontinuities arise 
in the flow for which the use of a Riemann solver can be an 
advantage. 

RH2D is a 2D mixed explicit/implicit second-order up- 
wind algorithm that also uses a staggered grid. It can treat 
radiation transport in the flux-limited diffusion approxima- 
tion, but in this paper we are only concerned with adiabatic 
discs. Its advection algorithm is based on the monotonic 
transport scheme by |van Leer| ( [19771 . 

FARGCQ is a 2D hydrodynamical code, using a polar 
grid centred on the star. It solves the Navier-Stokes and 
continuity equations, as well as the energy equation in a 
more recent version (Baruteau & Masset 2008a), which we 



use here. It is based on van Leer upwind algorithm, on a 
staggered mesh. Both FARGO and RH2D use the FARGO 
algorithm: in each ring i, at every time-step, the averaged 
azimuthal velocity v Vl i is computed. The ring is globally 
shifted by the corresponding number of cells for the consid- 



ered time-step length St : 



E{v Vt i-^-], where Sip is the 



elementary angle associated to a cell, and r the radius of 
the ring. Then, the advection is performed using the rem- 
nant azimuthal velocity in every cell v' v — v^ — rnr Sip I St. In 
rotating disks where \v' v \ <C v v , this enables a speed-up of 
the computation, and a lower numerical diffusivity becasue 
of the larger time step. 



4 ISOTHERMAL TYPE I MIGRATION 

In this section, we briefly review recent progress on Type I 
migration in the isothermal limit. This will prove helpful in 
understanding the general case, since very similar processes 
operate. 

One can linearise equations [I] and [2] and solve these 



rycansky & Pollack 19931. This yields a Lindblad torque 

(7) 



( Tanaka et al.||2002| their 2D result): 
r L /r = -3.2- 1.468a, 
witHU 



T = (q/h) £ p r p fip, 
and a corotation torque 
r c ,iin/r = 2.04- 1.36a. 



(8) 



(9) 



Note that the linear corotation torque is proportional to the 



radial gradient of specific vorticity (Goldreich & Tremaine 



19791, being zero for a = 3/2. It is also important to stress 
that these results were obtained with essentially no gravita- 
tional softening. 



It was shown in Paardekooper & Papaloizou (2009a I 



that corotation torques are non-linear in general, unless a 
very strong viscosity is applied. The linear corotation torque 
is replaced by non- linear horseshoe drag ( |Ward|[l991[ ): 



rHS/Po = 



3/3 \ 4 h 2 



(10) 



where a; s is the half-width of the horseshoe region, in units of 
r p . In the limit b — > 0, it was shown in [Paardekooper fc Pa-| 



paloizou (2009b I that x 2 = 1.68q/h, making the horseshoe 
drag 



Fhs/Fq 



3.18 -2.11a, 



(11) 



which is a factor of more than 3/2 larger than the linear 
corotation torque. The result that the non-linear corotation 
torque is larger than its linear counterpart also holds for 
non-zero gravitational softening. 

One can then combine equations [7] and [TT] to obtain a 
formula for the total torque: 

r/r = -0.02 - 3.578a, (12) 
which can be seen as a non-linear equivalent of the 2D for- 



mula of Tanaka et al. (2002): 



Tiin/ro 



-1.16 - 2.828a. 



(13) 



For a constant surface density disc (a = 0), inward migra- 
tion has slowed down by a factor of 100 through non-linear 
effects. However, these formulae are of limited use due to the 
lack of gravitational softening. In general b should be of the 
order of h to account for 3D averaging effects, which would 



lead to a different value of x a (Paardekooper & Papaloizou 



2009b I and a different Lindblad torque ( Paardekooper & 



Papaloizou 2009a I. The general conclusion that non-linear 



corotation torques can slow down Type I migration is still 
valid, however. 



5 A SIMPLE ADIABATIC MODEL 

In this section, we will construct a simple model describ- 
ing Type I migration in terms of a linear Lindblad torque 



http : //f argo . in2p3 . f r/ 



2 All torques presented in this paper will be normalised by Vg. 
Note that Tq is proportional to q 2 . 
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Figure 1. Total torque on a g = 1.26 ■ 10 — 5 planet embedded 
in an adiabatic disc (7 = 5/3) with a = 3/2 and j3 = 1, so that 
the corotation torque vanishes. Different curves denote different 
values of the softening parameter 6, and the dotted lines show the 
prediction of equation|14| Results were obtained with RODEO. 



plus the non-linear horseshoe drag. A first expression for the 
adiabatic horseshoe drag was proposed in Paardeko oper fc] 
|Papaloizou| ( |2008| , obtained by integrating the density per- 
turbation due to entropy conservation over the disc. In this 
approach, an assumption has to be made on the exact ge- 
ometry of the horseshoe region. Paardekooper & Papaloizou 
(20081 considered rectangular streamlines, and showed that 



the resulting torque is of the correct magnitude. However, in 
reality streamlines will not be rectangular, which can have 
a large impact on the torque. Here, we try to relax this as- 
sumption and take a different approach that allows us to 
combine the contributions of entropy and specific vorticity 
in a simplified way, and include the contribution of the Lind- 
blad torque. 



5.1 Linear torques 

Linearization of the two components of the equations of mo- 
tion, together with the continuity equation and the adiabatic 
condition (equations [TJ [2] and [3| yields a pair of first or- 
der ordinary differential equations (see Paardekooper &i Pa- 
paloizou 2008] their equations 14 and 15). We have solved 
these linear equations for different background surface den- 
sity and temperature profiles to obtain simple estimates of 
the Lindblad torque and linear corotation torque. 

Solving the linear equations (see [Paardekoo per fc Pa-| 
paloizou 2008) results in a linear Lindblad torque: 



7r L /r = -(2.5 + 1.7/3 -0.1a) 



b/h) 



(14) 



It was shown in Paardekooper & Papaloizou ( 2008 1 that the 



linear corotation torque in adiabatic discs is associated with 
singularities due to radial gradients in entropy and in the 




100 200 300 400 500 600 
t (orbits) 

Figure 2. Total torque on a q = 1.26 ■ 10 — 5 planet embedded in 
an adiabatic disc (7 = 5/3, h = 0.05) with different density and 
temperature profiles, for b/h = 0.4. Since the disc is adiabatic, 
the corotation torque saturates, leaving the Lindblad torque only. 
The dotted lines indicate the prediction of equation 1 1 4| Results 
were obtained with RODEO. 



quantity 
Ek 2 



(15) 



with k the epicyclic frequency, equal to Q in a Keplerian disc. 
The condition that the above quantity should be constant 
for the corotation torque to be zero is a generalisation of 
the condition that the gradient of specific vorticity should 
vanish, which applies in the strictly barotropic case. One is 
then lead to a two-term expression for the linear corotation 
torque; one proportional to £, and one proportional to 3/2 + 
(1 — 2/7)0 — 2/3/7. From our linear calculations, we found 



7rc,iin/r = 0.7 ( - + 



1 - 



2p 

1 



0.4 

7/ 7 J \b/h 

0.4 

b/h 

which can be written in terms of £ rather than /3 



+ 



2.2^ 



(16) 



7rc,iin/r = 0.7 ( - 



2* 



0.4 

b/h 



+ 



making the total linear torque 

Tim = Pl + Tcjin. 



(17) 



(18) 



We will compare this linear estimate to non-linear simula- 
tions at early times in section [5] 

For an isothermal disc (/3 = 0, 7 = 1, and therefore 
£ = 0), and b/h — 0.4, we have 



riao/ro 



-1.4 -0.6a, 



(19) 
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Figure 3. Lindblad torque on a q = 5 ■ 10 -6 planet embedded in an adiabatic disc (7 = 1.4) with h = 0.05. Solid lines indicate the 
prediction of equation |14| the dotted line indicates the 3D result of |Tanaka et al.| j2002| and the dashed line the 2D result from |Tanaka| 
|et al.| | [2002[ l. Symbols denote results from numerical simulations, obtained with FARGO. Left panel: fi = 0, for different values of a. Two 
values of b/h were considered: b/h = 0.6 (black solid line, black symbols), and b/h = 0.3 (grey solid line, grey symbols). Right panel: 
a = 1/2 with b/h = 0.6, for different values of 0. 



while a 3D calculation by Tanaka et al. ( 2002 1 resulted in 



,iso/ro 



-1.364 - 0.541a. 



(20) 



Therefore, our adopted value of the smoothing length gives 
reasonable agreement with fully 3D calculations in the 
isothermal limit. 

The dependence of the torque on softening can be quite 
complicated (see Paardekooper fc Papaloizou|2009a for the 
isothermal case). We have chosen for a simple power law 
scaling that is valid around b/h = 0.4, which is a reason- 
able value (see above). In Fig. [I] we show the total torque 
for a disc where the corotation torque vanishes, for different 
softening parameters. While equation |14| gives a good esti- 
mate for b/h = 0.6 (and b/h — 0.4, which is not shown), 
for more extreme values of b/h the simple estimate fails. 
However, these extreme values are not of interest physically, 
since b/h « 0.4 gives reasonable agreement with 3D results. 
For the sake of completeness, we note that the failure of 
equation [14] at small softening is not due to the failure of 
linearity, but due to the failure of the simple scaling with 
b/h of equation |14| 

Not only does the magnitude of the Lindblad torque 
depend on softening, also its dependence on a changes. This 
already can be appreciated by comparing the 3D isothermal 



result on the Lindblad torque from Tanaka et al. ( 2002 1 



3 /r = -2.34 + 0.099a, 



(21) 



to the 2D result given by equation [7] The coefficient of a 
changes sign between 2D (unsoftened) and 3D calculations. 
We also see this change when using smaller softening param- 



This is illustrated in Fig. [2] where we show the long-term 
evolution of the total torque for inviscid, adiabatic discs with 
various temperature and density profiles. Since there is no 
viscosity or heat diffusion, the corotation torque saturates, 
leaving only the Lindblad torque, which is then compared 
to equation |14| The agreement is very good for this value 
of b/h. Further experiments have shown that for smaller 
softening, the dependence on a is reversed, while for larger 
softening it is somewhat weaker. In all cases, however, the 
coefficient of a is small. 

This is illustrated in the left panel of Fig. [3] where we 
compare numerical results for two different values of b/h 
to equation |14| as well as to the formulae of Tanaka ct al. 
( |2002[ > . For b/h = 0.6, the coefficient of a is positive, and 
equation [14] gives good results for all values of 01. This also 
holds for b/h = 0.4. For smaller values of b/h, the trend 
with a reverses, and equation[l4]deviates from the numerical 
result by approximately 15% for a — 0. 

In the right panel of Fig. [3] we show the trend of the 
Lindblad torque with 0. We have also reconstructed the /3- 



dependence of the linear results from Tanaka et al. ( 2002 1 



(their tables 1 and 2). Their 3D result gives a temperature 
dependence that is less steep than we find from our 2D sim- 
ulations. The numerical results are in very good agreement 
with equation |14| 



5.2 Horseshoe drag 

In this section, we present a simple model for the non-linear 
corotation torque, the horseshoe drag, in the presence of 



eters. Equation 14 gives good results for b/h ~ 0.4, however. 



entropy and vortensity gradients. Following Ward ( 1991 1 we 



6 S.-J. Paardekooper, C. Baruteau, A. Crida and W. Kley 



A A A A 



3.20 



3.15 



3.10 



3.05 



WW V V V V 



A M A A AAA 



V V \/ \/ 



0.96 0.98 



1.00 1.02 1.04 

r 



Figure 4. Schematic overview of streamlines near the horseshoe 
region, with the planet indicated by the black dot. The separatrix 
is colored black, and the thick curve indicates the location of the 
entropy discontinuity. The top and bottom of the figure can be 
thought of as T2 and IT , respectively. This picture applies before 
any material that has made the turn comes back on the other side 
of the planet. 



consider the torque produced by material on streamlines un- 
dergoing horseshoe turns. We consider a region TZ interior to 
the two separatrices, separating the horseshoe region from 
the rest of the disc, and bounded by two lines of constant 
if, 1\ and T2 on the trailing and leading sides of the proto- 
planet respectively (see |Paardekooper fc Pa paloizou 2009a, 
and also Fig. j4j). These boundaries are supposed to be suffi- 
ciently far from the protoplanet that the corotation torque 
is determined within. Assuming a steady state, this torque 
may be obtained by considering the conservation of angular 
momentum within TZ written in the form 



IL 



E ) r,!,s,lr 



1 r 2 



Fdr 



(22) 



with F = E(j — j p )(fi — Q p )r (see Paardekooper & Pa- 
paloizou 2009a). Here j = rv v is the specific angular mo- 



mentum and j p is j evaluated at the orbital radius of the 
protoplanet. Assuming symmetric horseshoe turns (for a val- 



idation see section 5.2.2 below), we have 



2r- p f 
Jo 



(F - F )dx, 



(23) 



where Fo equals F in the unperturbed disc and x — (r — 
r p )/r p . Assuming Keplerian rotation, and using first order 
expansions for j — j p and Q — fl p , we can approximate 



F-Fo 



-r 3 T Q 2 t 2 

4 Ijq 



E 



(24) 



Below, we discuss some physical arguments that allow us 
to relate the state (density, pressure and velocity) after the 
turn to the initial state. 



5.2.1 Pressure equilibrium 

First of all, it is important to note that the disc will always 
try to maintain pressure balance: material that has executed 
a horseshoe turn should still be in pressure equilibrium with 
its surroundings. This assumption was not made in the orig- 
inal barotropic model ( |Ward|1991[ ), where the density (and 
therefore the pressure) was allowed to change while keeping 
the rotation profile fixed, similar as in equation |24| In real- 
ity, the disc will change its rotation profile in order to retain 
pressure balance. However, this adjustment after the turn 
does of course not affect the torque. We can therefore use 
equation |24[ where only density changes are considered, to 
obtain the torque on the planet, but we have to keep in mind 
that the actual state of the fluid after the turn may well be 
different. Below, we work out the density changes due to en- 
tropy conservation and vortensity evolution. While entropy 
conservation can work directly on the density, because the 
pressure can be kept constant, the vortensity evolution will 
mainly affect the rotation profile, just as in the barotropic 
case. To obtain the torque, however, it is again sufficient to 
consider changes in surface density only, and use equation 



5.2.2 Mass conservation 

One important constraint that has not been discussed so 
far is conservation of mass, which states that the amount of 
mass that goes into the horseshoe turn at x < must come 
out of the turn at x > 0. If material between x = xo < and 
x — will make the turn and come out of the turn between 
x — and x = x\ > 0, then we must have that 



/' 



S(n - fl p )rdr = 0. 



(25) 



It is clear that the assumption of pressure balance, together 
with mass conservation, will lead to an asymmetry in the 
horseshoe leg (i.e. xi 7^ — xo). Since the effect works in the 
opposite direction for the other horseshoe leg, the end result 
is that one leg will appear wider than the other. Although 
this in general can affect the torque, we show below that the 
impact on the torque is usually quite small. 

For simplicity, we take the gradients of entropy and spe- 
cific vorticity to be zero. There is still an asymmetry in this 
case due to curvature, and this will give us an estimate of 
the importance of this effect. We then have E = £0 and 
Q — Qk (ignoring the radial pressure gradient), and using a 
first-order Taylor expansion of — Q p and Eq, we have 



(1 — ax + x)xdx — 0. 



(26) 



Writing x± = —xq + 8x, and keeping only terms that are 
first order in Sx, we can solve for Sx: 



5x= - (a 



1) l n 



(27) 



Since Sx <C xo (because xo 1), the combined effect of 
both horseshoe legs on the torque, which scales 
change of a factor (1 + 2Sx/xq), which can be of the order of 
10% for considerable gradients in density. When gradients 
of specific vorticity and entropy are present, they can also 
contribute to the asymmetry (for example through equation 
28 below). For gradients in entropy and vortensity that are 
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Figure 5. Top panel: source term in the vortensity equation, 
in units of f2p/£ p . Bottom panel: relative change in vortensity 
(u;/£ — ujo /So) / (u>o /So) • The dashed curves show the location 
of the separatrix. Both panels show results for a q = 1.26 ■ 10 _o 
planet, embedded in an h = 0.05 disc with a = 3/2 and j3 = —3/2 
after 10 orbits. A large value for softening parameter was used, 
b/h = 2, to reduce the influence of the wakes on the coorbital 
region. Results were obtained with RODEO. 



not too large, this effect is small. The assumption made in 
equation [23] above that the horseshoe turns are symmetric 
should therefore give reasonable results. 



5.2.3 Entropy conservation 

In pressure equilibrium, changes in entropy are directly re- 
lated to changes in density. Consider material that has made 
a horseshoe turn from —x to x (we assume symmetry at 
this point; see section 5.2.21. The old state at x is given 
by So, po and so, the new state by S, p, s. Pressure bal- 
ance dictates that p = po, while entropy conservation gives 
s = so(l + 2£a;), which makes the density S = (po/s) 1 / 7 , or: 



S = 



So (l 

S 



2«z 



< x < x s 
otherwise. 



(28) 



5.2.4 Specific vorticity 

In a barotropic disc, specific vorticity is conserved along 
streamlines. We can, analogous to the entropy case discussed 
above, write down an expression for the change in specific 
vorticity after the turn: 

oj/S — wo/So _ 2 ^dlog(o;o/So) 



(29) 



dlogr 

and zero otherwise. Here, us — V x v denotes 



ui / S 

for < x < X., 
the vorticity. 

In a non-barotropic disc, specific vorticity, or vortensity, 
is no longer conserved along streamlines: 

D / ui \ _ VS x Vp _ Vs x Vp 
Dt VS/ ~ 



S3 



7S 2 



(30) 



In a barotropic disc, in which p = p(S), VS and Vp are 
parallel everywhere, which makes the source term in above 
equation zero, with the result that vortensity is conserved 
along streamlines. Note, that for a non-barotropic disc, in 
regions where the dens ity and pressure are smooth, the right- 
hand side of equation 30 is small, since p ~ /i 2 S with h <C 
1. However, at the outgoing separatrix, material that has 
made a horseshoe turn meets disc material that still has the 
unperturbed value of the entropy. At this specific streamline, 
a large entropy gradient exists (it is formally infinite across 
the separatrix) perpendicular to the flow (see Fig. The 
change in vortensity induced by the associated source term 
gives rise to an additional entropy-related torque. 

This is further illustrated in Fig. [5] where we show the 
2D distribution of the specific vorticity (bottom panel), to- 
gether with the source term in equation [30] (top panel). The 
background surface density profile is such that the vorten- 
sity is constant, initially. All structure seen in the bottom 
panel of Fig. [5] is therefore due to the source term depicted 
in the top panel. It is clear that this source term only acts on 
the outgoing separatrix, where advection of entropy gener- 
ates an entropy discontinuity. The pressure gradient entering 
equation [30] is due to the hydrostatic envelope of the planet, 
hence the source term is localized around the planet. En- 
tropy advection along the horseshoe bend does not change 



A similar equation holds for the other horseshoe leg. 



the pressure ( Paardekooper & Papaloizou 2008 Baruteau 
fc Masset|2008a 1, so the only other possible pressure gradi- 
ents are due to the wakes (which play no role for low-mass 
planets, where x s < h), or due to a global radial pressure 
gradient. The latter can indeed contribute to the vortensity 
source, but it is easy to see that the effect will be symmetric 
in both horseshoe legs, and therefore this does not result in 
a torque onto the planet. In Fig. |5j the initial pressure was 
taken to be constant (/3 = —a = —3/2) for clarity. 

It can be seen from Fig. [5] that the lines of constant 
specific vorticity are not exactly parallel to the separatrix, 
moving slightly away from the planet's orbit near <p = 2.0 
and tp — 4.5. This is due to the adjustment of the disc 
to maintain pressure equilibrium. The propagation of the 
vortensity discontinuity triggers a pressure wave that can be 
clearly identified especially at early times. It can also be ob- 
served in barotropic discs whenever there is an initial radial 
gradient in specific vorticity. The presence of the vortensity 
discontinuity makes this pressure wave more apparent. 

From now on, we assume that vortensity is conserved 
everywhere except along the outgoing separatrix. Consider 
the integration of equation |30] along the outermost stream- 
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line of the horseshoe region. The only entropy gradient of 
importance is perpendicular to the streamline after it has 
encountered the planet (see above). We therefore write equa- 
tion [30] as: 

D /u\ VxsV||P Vj_flV||II 

Dt VE/ 7 E 2 s 7 Es ' ( ' 

where Vi and Vm indicate the component of the gradient 
perpendicular and parallel to the streamline, respectively, 
and II is the fluid enthalpy. The orientation is such that we 
take the gradient of IT in the direction away from the planet, 
which means that we have to take the gradient of s into the 
horseshoe region. 

We take the gradient of s into the horseshoe region to 
be infinite, formally: 

V_ls/s = 2£x s 6(x-x B ). (32) 

It is easy to see that this models a jump in entropy across 
the separatrix that has the correct magnitude. 

Noting that factors involving s are constant along the 
streamline, and assuming E m E p and taking the velocity 
along the streamline to be v = vr p Q p x s , with v a constant, 
we can integrate equation [31] from the turn at x = near 
the stagnation point to a point far away from the planet to 
end up with: 

2£ ritum - n (a; B 



a (e) 



-S(x- 



(33) 



where n tur n is the enthalpy at the location of the turn, 
near the stagnation point. Along this streamline, the veloc- 
ity varies from at the stagnation point to v = 3x s r p Q, p /2 
far away from the planet; therefore we need < v < 3/2. 
The choice of v basically depends on the exact geometry 



of the horseshoe region (see below). Equation 33 then gives 
the the vortensity production at the outgoing separatrices, 
a process that does not operate in barotropic discs. 

Note that in the case of multiple stagnation points 



close to the planet (IMasset et al. 120061 IPaardekooper & Pa- 



paloizou 2009b see also Fig. 131, there is always a single 



stagnation point where the entropy discontinuity starts for 
both horseshoe legs (the stagnation points below the planet 
in the top panels of Fig. 131. We therefore do not have to 



make any additional assumptions on the detailed flow topol- 
ogy close to the planet. 

To make further progress, we now assume that the pres- 
sure (or enthalpy) structure is not significantly changed from 
the barotropic case, or, equivalently, the case with £ = 0. 
For £ ^ 0, advection of entropy leads to changes in en- 
thalpy, which is then discontinuous across the separatrix. It 
is difficult to see which value of II to take (inside or out- 
side the separatrix) in that case. Note, however, that the 
jump in IT is of order x B and therefore small. It should not 
affect the pressure structure near the outgoing separatrix. 
This has been verified using numerical simulations. We can 
then find fl( urn — FIq by using the Bernoulli invariant E (sec 



) 2 + n + $ p - lr 2 p Q. 2 p x\ 



(34) 



Paardekooper fc Papaloizou|2009b| | : 
1 

2' 

Considering two points on the same streamline, one at the 
turn (where x = and Q, = fi p ) and one far away from the 
planet (where x = x s and $ p =0), we have: 

+ $ P ,turn = -\x 2 s r p Vt 2 p + n (a; s ). (35) 



n, 



8 



Using equation |35| in equation |33| we end up with 



A(|)W 



2^ ^ 

«7 Ep 



x. 

d 8 6 



S(x — x s ), 



(36) 



with d — \/jrturn — r p| 2 / r P + b 2 . Numerical simulations for 



isothermal discs (Masset et al. 2006 Paardekooper & Pa- 



|paloizo u 2009b) show that the stagnation point is located 
roughly 3/2 softening lengths away from the planet, mak- 
ing d = J 13/46. In adiabatic discs, the situation is slightly 
different. As was shown in Masset et al. (20061, the width 



of the horseshoe region is directly related to the perturbed 
value of the Bernoulli invariant at the stagnation point 
E = II turi) + $ p ,turn: 



r p Q p 



-E' 



(37) 



On the other ha nd, it was shown in [Paardekoo per fc Pa-| 



paloizou (2009b I that a; s oc I/7. Therefore, there should 



exist a direct relationship between the location of the stag- 
nation point and 7. We should have E' oc 1/^/7, and, as- 
suming for simplicity that both the perturbed enthalpy and 
the planet potential at the stagnation point have the same 
dependence on 7, this leads to 



r - r t „ 



= \fnrj - lb, 



(38) 



for some constant n. It was noted in Paardeko oper fc Pa-| 
paloizou ( 2009b I that the location of the stagnation point is 
essentially determined by the Lindblad wakes, which makes 
it very difficult to model. The best we can do is fix the con- 
stant n so that we recover the isothermal result for 7=1. 
This means we take n = 13/4, and we will see in section 6.3 



that this gives reasonably good results for 7 > 1. We then 
have 



d- 



(39) 



Away from the outgoing separatrices, the specific vor- 
ticity source term is small, and we can assume conservation 
of vortensity. We then have for the total vortensity: 



E 



So 



d log r 



< x < x s 
otherwise. 



(40) 



The term in parenthesis results from conservation of specific 
vorticity, that is also present in barotropic discs. 

5.2.5 Total horseshoe drag 



We now add the contributions of entropy (equation 28 \ and 



vortensity (equation 40 1 to the density perturbation, which 
gives 



E 



E 



-2-x + 2 



-A 



Wo 



(41) 



The vortensity-related perturbation (the second term on the 
right-hand side) is the same as in the barotropic case. The 
entropy-related density perturbation (the first and the last 
term on the right-hand side) is caused by density structures 
produced by material conserving its entropy, bound to the 
horseshoe region, plus an additional component linked to 
the production of vortensity at the outgoing separatrices. 
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We stress again that the contribution of the vortensity to 
the density perturbations (the last two terms in the equation 
above) will affect the rotation profile rather than the density. 
However, the torque exerted on the planet will be the same, 
just as in the barotropic case. 

In order to find the torque, we can now perform the 
integral in equation |23[ using equations |24| and |4 1 1 yielding 



4q2 4 



a + 



vdx 2 



+ 1 



(42) 



For £ = 0, we recover the barotropic result of |Ward| 
(19911. Note that although x s is well-defined for b — > 



\ Paardckoopcr & Papaloizou 2009b), the contribution of the 
entropy discontinuity diverges if d — > as well. 

If we write x s = Cx/g/VV 74 , with C = C(b/h) 



I jPaardekooper fc Papaloizou|2009b[ ), then we have: 



jT c , hs /r = \c 4 



>y\vC 2 d 



+ 1- 



Note that the component due to a radial entropy gradient 
can easily overpower the contribution from the vortensity 
gradient. Note also that for fixed b/h (which makes C a 
constant as well, as long as a; s < h), the horseshoe drag 
r Cj hs scales as q 2 /h 2 , just as the linear torque. 

Numerical simulations and analytical arguments indi- 
cate that the geometry of the horseshoe region is the same 
for all low- mass planets, as long as x s < h (Masset et al 



2006 Paardekooper & Papaloizou 2009b I. Therefore, one 



7r c ,hs/ro = 1.1 



0.4 

b/h 



+ 



g 0.4 

7 b/h 



10.1 




2.2 



-10 



(43) 
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choice of v should suffice. We have obtained good agree- 
ment with numerical simulations using v — 1.0. We have 
measured the horseshoe width x s to be 



(44) 



i.e. C — 1.1 for 6 = OAh. We note that the scaling with b/h 
breaks down for small softening (b/h < 0.3). Equation |43| is 
then completely determined: 



Figure 6. Total torque on a q = 1.26 ■ 10~ 5 planet (with 
b/h = 0.4) embedded in an adiabatic disc (7 = 5/3, h = 0.05) 
with different density and temperature profiles. The dotted lines 
indicate the prediction of equation [l8] (middle two lines) and the 
result of equatio r|47| (top and bottom lines). Results were obtained 
using RODEO. 



middle panel of Fig. [7] for 7 = 5/3. Positive torques, and 
therefore outward migration, are readily obtained for f3 > 0, 
i.e. for temperature profiles that decrease outward. A glance 
at equation [14] reveals that this is completely due to the 
entropy-related corotation torque, since the Lindblad torque 
becomes more negative with increasing /3. In the top panel of 
Fig.[7|we show equation[47]for 7 = 1.4, a value often adopted 
for protoplanetary discs (Paardek ooper fc Papaloizou|200"8 



5.3 Total torque 

The total torque in the non-linear regime (t >~ 2 orbits 
in an inviscid disc), before saturation sets in, is given by 
( Paardek ooper fc Papaloizou|2009a I 

r = r L + r c , hs , (46) 

with Pl given by equation |14| and r Cj hs given by equation 



We now take b/h — 0.4, which makes the total torque: 
7 r/r = -2.5 - 1.7/3 + 0.1a + 1.1 f | - aj + 7.9^, (47) 

where the last two terms describe the non-linear corotation 
torque. This equation is compared to numerical simulations 
in Fig. [6j showing remarkably good agreement between our 
simple model and fully non-linear hydrodynamical simula- 
tions. We will defer a detailed numerical analysis to Sect. 

El 

Equation [47] is shown as a function of a and j3 in the 



Kley fc CridapOOS) ). For this lower value of 7, the entropy 
gradient depends stronger on the temperature gradient (see 
(45) 

equation |5|, resulting in a steeper slope for the zero-torque 
line. Since protoplanetary discs are expected to have f3 > 
in most parts, this indicates that outward migration is a 
serious possibility. 

In the bottom panel of Fig. [7] we return to 7 = 5/3, 
but use a larger value of C, C = 1.3, that would correspond 
to higher mass planets that are able to push against the 
Lindblad wake to take the stagnation point close to (r, tp) — 
(r p , 7r) ( Masset et al.|20"06 l. Formally, one should use d = b 
as well, but we have found that the stagnation point never 
actually reaches the planet. Effectively, one should use C < 
1.3 in combination with b < d < \/5fc, but C = 1.3 in 
combination with d = \* r 5b gives reasonably good results 
(see Sect. |6.6~|. 



5.4 Locally isothermal limit 

It is straightforward to obtain the isothermal limit of equa- 
tion 



47 



by setting [3 = and 7 = 1 (and therefore £ = 
0), which then leaves the linear Lindblad torque plus the 
vortensity-related horseshoe drag: 



r iso /ro = -0.85 - a, 



(48) 
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Figure 7. Total torque, as given by equation |47| in units of To/7. 
The solid line indicates where the total torque is zero. The black 
dot indicates the Minimum Mass Solar Nebula, having /3 = 1 and 
a = 3/2. Top panel: 7 = 1.4, middle panel: 7 = 5/3. The bottom 
panel is the same as the middle panel, but for C = 1.3, valid for 
higher mass planets that obtain the maximum value of x B . 




Figure 8. Total torque on a q = 5 • 10~ 6 planet (with b/h = 0.6) 
embedded in a locally isothermal disc with h = 0.05. Diamonds 
denote numerical results obtained with FARGO, squares the re- 
sult of equation |49| Corresponding results are connected with a 
line. For the top row (7 runs), a = 1/2, for the middle row (4 
runs) a = 3/2, and for the bottom row (2 runs) a = 5/2. 



where we have used b/h = 0.4. Compared to the 3D linear 



result (equation 20 1, migration has slowed down by approx- 
imately 50% for a = 0, but the direction is inward for all 
realistic surface density profiles. 

A different approximation that is often used is the lo- 
cally isothermal limit, which means solving the isothermal 
equations but with a radially varying sound speed (or, equiv- 
alently, temperature). One arrives at this limit by taking 
•y 1, and, crucially, invoke infinitely efficient thermal dif- 
fusion. This effectively takes the entropy-related horseshoe 
drag into the linear regime. The total torque then consists of 
the linear Lindblad torque, the linear entropy-related corota- 
tion torque plus the non-linear vortensity-related horseshoe 
drag: 



rioci S o/ro = -(2.5 - 0.5/3 - 0.1a) 



1.4(3 



0A\ 
b/h) 



1.1 



b/h) 

(-) 
{b/h) 



(49) 



Note that the temperature dependences of the Lindblad and 
corotation torque work against each other, leaving a total 
torque that is less sensitive to temperature variations. 

In Fig. [8j we compare equation [49] to numerical results 
obtained with FARGO, showing good agreement over a wide 
range of (3. Note that, contrary to the adiabatic case, a neg- 
ative temperature gradient works in favour of inward mi- 
gration. This is due to the relatively weak dependence of 
the linear corotation torque on /3, compared to the Lindblad 
torque. We have found inward migration for all reasonable 
values of a and (3. 

However, the good agreement as seen in Fig. [8] is not 
the whole story, as can be seen from Fig. [9] While the solid 
curve denotes a case of constant specific vorticity, a non- 
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Figure 9. Total torque on a q = 1.26 • 10 -5 planet (with b/h = 
0.4) embedded in a locally isothermal disc (7 = 1, h = 0.05, 
P ^ 0) with different density and temperature profiles. The dotted 
lines indicate the result of equation |18| Results were obtained with 
RODEO 



linear rise in the torque can be observed. This is due to the 
source term in the vorticity equation. It is important to note 
that the analysis presented above for adiabatic discs is not 
valid in the locally isothermal case, since entropy (which 
would correspond to in this case) is not conserved along 
a streamline. One would get a source term proportional to 
dp /dip d,Cs/dr, the effect of which will strongly depend on 
the geometry of the horseshoe region. Note, however, that 
the impact of the source term is quite small even for a steep 
temperature gradient /3 — 2. We note that the effect is neg- 
ligible for |/3 1 < 1, and for /3 = 2 comparable to the non- 
linear effect associated with barotropic horseshoe drag (see 
the dashed curve in Fig.[9j|. A locally isothermal disc is nev- 
ertheless an interesting case, since it would correspond to a 
part of the disc that can cool very efficiently, i.e. the outer 
parts of protoplanetary discs. 



6 NUMERICAL RESULTS 

In this section, we will try to validate the simple model by 
means of numerical simulations. 
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Figure 10. Total torque on a q = 1.26 ■ 10 planet (with b/h = 
0.4) embedded in an adiabatic disc with 7 = 5/3, h = 0.05, 
o = and /3 = 1. Different line styles denote different codes, 
while different line colors indicate whether the planet is kept in 
the middle of a grid cell or on the edge. The dotted line indicates 
the result of equation |47| The resolution used is Ar/r p = Aip = 
0.0025. 
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6.1 Code comparison 

We first check whether our independent numerical methods 
give similar results. In Fig. |10| we compare results for the 
three different numerical methods on an adiabatic disc with 
7 = 5/3, h — 0.05, a = and /3 = 1. The grey curves 
denote models for which the grid is chosen so that r/r p = 1 
lies on the edge of a grid cell. For these models, all codes 
agree nicely on the final torque, with each other as well as 
with equation |47| They also agree on the linear part of the 
torque (for t < 2). RODEO shows a slightly faster rise in 



Figure 11. Total torque on a q = 5T0 -6 planet (with b/h = 0.6) 
embedded in an adiabatic disc with h = 0.05 and 7 = 1.4. Dia- 
monds denote numerical results obtained with FARGO, squares 
the result of equation |46| Corresponding results are connected 
with a line. For f = —0.7 and £ = 0.8, three values of a were 
considered, a = {0.5, 1.5,2.5}, from top to bottom. For £ = —0.4 
and £ = 0.4, we show results for a = 0.5 and 1.5 (from top to bot- 
tom), while for all other values of £ results are shown for a = 0.5 
only. In all cases the temperature gradient can be found from 
P = £ + 0.4a. 
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Figure 12. Same as Fig. |10| for different radial resolutions. All 
results were obtained with FARGO. 



the torque. Since this difference arises as soon as non-linear 
effects set in, it probably originates close to the separatrix, 
where the vorticity source term plays a major role. 

In Fig. we compare the difference between equation 
[46] and numerical results obtained with FARGO for different 
entropy gradients. The relative difference is well within 20% 
except for the cases where £ = 0.8. Results obtained with 
RODEO show similar good agreement. Therefore, not only 
do all codes agree on the torque, they also agree very well 
with our simple analytic model over a large range of a and 
(3 within 30%. 

Results displayed in Fig. [10] were obtained at a reso- 
lution of Ar/r p = A(y3 = 0.0025. For C = 1.1, we have 
x s = 0.015, so the half-width of the horseshoe region is re- 
solved by 6 grid cells. Lowering this to 4 cells had little effect 
on the torques, suggesting that they are converged at this 
resolution. There is however the interesting difference be- 
tween the black and grey curves in Fig. |10| which we discuss 
next. 



6.2 Planet position 

The only way in which the simulations indicated by the 
grey and the black curves in Fig. [lO] differ is the position 
of the planet on the grid. While for the results obtained 
with RODEO, this has almost no effect on the torque, both 
FARGO and RH2D show a difference of approximately 50% 
between the two planet positions. This is a numerical effect 
that has to do with the sharp gradient in specific vorticity 
that arises at the outgoing separatrix under influence of the 
entropy-related source term. Such a sharp feature is difficult 
to handle numerically. It also appears in locally isothermal 
simulations, but is absent in fully isothermal runs, where 
there is no source term for the specific vorticity. The use of a 
Riemann solver, which is specifically designed to handle dis- 
continuities, largely evades this problem. Being a numerical 
artefact, it also disappears at higher resolution, which is il- 



lustrated in Fig. 1 12| After doubling the radial resolution, the 
grey curve becomes indistinguishable from the black curve. 
The torque then does not show the slow decay after 10 or- 
bits, and settles close to the value that is obtained by putting 
the planet on the edge of a grid cell at low resolution. 

The amplitude of this numerical effect strongly depends 
on the detailed geometry of the horseshoe region close to 
the stagnation point, which in turn depends on background 
gradients of temperature and density. We observed that it 
is virtually absent for a disc with a = 3/2 and /3 = —3/2. 



6.3 Streamline analysis 

The exact geometry of the horseshoe region is of signifi- 
cant importance in determining the corotation torque. In 
isothermal discs through x s only, but for adiabatic discs also 
through the location of the stagnation point r sta g. There 
is a direct relation between x B and r ata g (|Masset et al 



20061 IPaardekooper & Papaloizou||2009b|). In Fig. 113 



show the streamlines close to the planet for isothermal mod- 
els (left panels) and adiabatic models (right panels) for a 
q = 1.25 ■ 10~ 5 planet (top panels) and a q = 1.08 • 10~ 4 
planet (bottom panels). The top left panel is consistent with 
the findings in Masset et al. ( |2006[ ), with the stagnation 
point located approximately 36r p /2 from the planet. For an 
adiabatic model with the same disc parameters, the stagna- 
tion point shifts further away from the planet, to approxi- 
mately 2br p . It was checked that an adiabatic model with 
7 = 1.01 looks similar to the isothermal model. We conclude 
that the simple scaling of equation [39] gives a reasonable es- 
timate for the location of the stagnation point for different 
values of 7. We note again that it is very difficult to model 
the position of the stagnation point in a simple way, be- 
cause it depends on the way the wake is able to influence 
the corotation region ( Paardekooper fc Papaloizou|200"9b \ . 



6.4 Softening 

In this section we discuss models that use different gravi- 
tational softening parameters. Although we have argued in 
that using b/h 



Sect. 



5.1 



0.4 gives linear torques that are 
in agreement with 3D linear theory, in is not clear at present 
what value of b/h will reproduce the 3D non linear torque 
(or horseshoe drag). It is therefore important to investigate 
the behaviour of the torque as a function of b/h. 

For isothermal models, it has already been noted that 
there is a strong dependence on b/h, due to the C 4 scal- 
ing of the vortensity-related horseshoe drag (Paardekooper 
& Papaloizou 2009a). Smaller values of b/h give stronger 
corotation torques, but since x s is finite in the limit b — > 
( jPaardekooper &: Pa paloizou 2009b| the horseshoe drag is 
well-defined. The first term of the entropy-related horseshoe 
drag (see equation 43 1 is proportional to C 2 h/b, and there- 
fore formally diverges for b — > 0. In our simple model, this 
stems from the stagnation point being located where the 
planet potential diverges. In practice, there will always be a 
finite distance from the stagnation point to the planet, and 
furthermore the planet potential should not diverge. For ap- 
propriate values of d and therefore <3>turn equation |43| should 
still hold in the limit b ^ 0. This limit is impossible to reach 
through numerical simulations, of course, since one always 
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Figure 13. Streamlines close to the planet embedded in a disc with a = and /3 = 0, for 7 = 1 (isothermal models, left panels) and 
7 = 5/3, right panels. Top panels: q = 1.26 ■ 10 — 5 , bottom panels: q = 1.08 ■ 10 -4 . Solid curves denote streamlines belonging to the 
horseshoe region, while dashed curves denote streamlines from the inner and outer disc. Results were obtained with RODEO. 



must resolve the gravitational potential in hydrodynamic 
simulations. 

We have varied b/h between 0.2 and 0.6, and the re- 
sults are displayed in Fig. |14| Black lines indicate isothermal 
simulations for a — 0, which show a less negative torque 
for smaller softening parameters. Therefore, the horseshoe 
drag, being positive and stronger for smaller values of b/h, 
is more than able to compensate for the Lindblad torque, 
which is more negative for smaller b/h. This effect, due to C 
being larger for smaller values of b/h, is nothing compared 
to the softening dependence of the entropy torque (see the 
grey curves in Fig. |14[ |. The total torque varies by an order 
of magnitude in this range of b/h, and this is including the 



more negative Lindblad torques for smaller b/h. Clearly, b/h 
is a crucial parameter for the entropy-related torque. 



The dotted lines in Fig. 
equations similar to equation 



have been obtained from 
for appropriate values of 



C and d. We have found that a constant value for d/b gives 
good results, i.e. the distance between the stagnation point 
and the planet is a fixed number of softening lengths, ir- 
respective of the value of b/h. The scaling of equation |44| 
breaks down, however, for b/h < 0.3, and we have measured 
C = 1.26 for b/h = 0.2. This gives a good match to the sim- 
ulation for the total torque. We note that in principle, one 
could come up with a more complicated functional form for 
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Figure 14. Total torque on a q = 1.26-10" 5 planet embedded in a 
h = 0.05 disc with a = 0, for different values of b/h, for 7 = 1 and 
(3 = (black curves), and for 7 = 5/3 and (3 = 1 (grey curves). 
Dotted lines indicate the result of equation |47| for appropriate 
values of C and d. Results were obtained with RODEO. 
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Figure 16. Perturbed vorticity (w — ljo)/u>o due to the presence 
of a q = 1.26-10 planet n a disc with h = 0.05, a = and /3 = 1 
after 20 orbits. The strong entropy gradient at the outgoing sep- 
aratrix leads to the appearance of an anti-cyclonic vortex, visible 
as the circular shaped vorticity minimum at (r, ip) = (0.985, 5.2). 
Results were obtained with RODEO. 
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Figure 15. Total torque on a q = 1.26- 10 — 5 planet embedded in 
a h = 0.05 disc (7 = 5/3) with a = 0, /9 = 1, for different values 
of b/h. Results are shown for RODEO (squares) and FARGO 
(diamonds) . The result of equation |46| is indicated by the solid 



C(b/h), using the results from Paardekooper & Papaloizou 
( |2009b| ). 

Although our numerical methods agree on the strong 
dependence of the torque on b, for small softening the dif- 
ferences become larger. This is illustrated in Fig. |15[ where 
we compare the total torque for different softening parame- 
ters for RODEO and FARGO. For b/h > 0.4 the agreement 



is very good, while for b/h — 0.2 the difference is approx- 
imately 30%. Note that the source term in the vortensity 
equation becomes very strong at small softening, leading to 
a strong contribution from the entropy discontinuity at the 
outgoing separatrix. This is a very challenging situation for 
numerical methods, and it is not surprising that the differ- 
ences between the methods are larger in this regime. 

The solid curve in Fig. [15] denotes equation |46| with- 
out the correction applied in Fig. [14] for the breakdown of 
the scaling of C with b/h. Therefore, it predicts too large 
a torque at small softening (about 20%). Overall, however, 
the agreement is very good. 



6.5 Vortex formation 

Vortices in protoplanetary discs can form as a result of the 
Rossby wave instability ( Lovelace et al.|1999 1 . This instabil- 
ity is usually discussed in the context of giant, gap opening 
planets ( |Li et al.|2005||de Val-Borro et al.|2007[ ), for which it 
was observed in most numerical codes for inviscid discs (de 
I Val-Borro et al.||2006| |. Although the low-mass planets dis- 
cussed in this paper are not massive enough to significantly 
perturb the surface density, the strong entropy gradients 
that exist at the outgoing separatrix can lead to vortex for- 
mation. This was observed in almost all simulations with a 



significant initial entropy gradient (see Fig. 16 for a typical 
example). Baruteau & Massct (20 08a[ ) also reported the ap- 
pearance of a vortex in the same context. When the source 
term in the vorticity equation acts to decrease the vorticity, 
which is always true for one of the horseshoe legs, there is 
the possibility of forming an anti-cyclonic vortex. Since the 
formation of this vortex occurs after the turn, it does not 
affect the torque. This may not be the case when the vortex 
interacts with the planet when it reaches the opposite side; 
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Figure 17. Total torque on planets of different mass (using b/h = 
0.4), embedded in an adiabatic 7 = 5/3 disc with h = 0.05, a = 
and P = 1. The dotted lines denote equation [47] for C = 1.3 (top) 
and C = 1.1 (bottom). Results were obtained with RODEO. 



it may then affect the partial saturation of the corotation 
torque. This will be discussed in a forthcoming work. Note 
that a small kinematic viscosity is enough to kill the vortex 
before it reaches the opposite side of the planet. 



6.6 Higher- mass planets 

Although the main focus of this paper in on l ow-mass plan- 
ets, for which a; a < h, or, equivalently, q < h 3 ( Paardekooper 
|fc Pa paloizou 2009b), there exist an interesting intermedi- 
ate class of planets with masses of approximately 20 - 50 
Mm that experience a boost in the corotation torque ( Mas- 



set et al.|2006| . It was shown in Paardekooper & Papaloizou 
( 2009b I that this is due to the fact that for x s > h, Lind- 



blad torques are less effective in affecting the shape of the 
horseshoe region. As a result, the width of the horseshoe re- 
gion increases, and the stagnation points move closer to the 
planet. Both effects will enhance the horseshoe drag, in the 
isothermal case ( |Masset et al.|[2006[ ), but especially in the 
adiabatic case. 

In Fig. |17[ we show the total torque on planets of dif- 
ferent mass. Note that since To oc q 2 , curves for planets 
for which x s < h would fall on top of each other when the 



x-axis is rescaled in an appropriate way (Paardekooper & 
|Papaloizou|2009"a| ) . This is expected from linear theory, but 
this also holds when the non-linear horseshoe drag is incor- 



porated, because x a o c ~Jq ( Masset et al 
& Papaloizou||2009b|). We 



2006[ 
found that for qjH 



Paardekooper 



> 0.2, the 
(2006J sets in. This can 
be seen in Fig. [17] from the solid and dashed curves. For 
q = 10~~ 4 = O.S/i^the maximum torque is reached. Beyond 
this mass, gap formation starts to play a role, lowering the 
mass and the opacity of the corotation region, which reduces 
the effect of the corotation torque. This was also observed 



The maximum torque in Fig. 17 achieved for q — 10~ , 
is in good agreement with equation |47| but with C = 1.3. 
As mentioned in Sect. |5.3[ since the stagnation point does 
not actually reach the location of the planet (see the bot- 
tom panels of Fig. 



13) 



47 



The overall result is 
1.3 (see Fig. 17 1, which is good enough 



C < 1.3, but this will be at 
least partly compensated by the d ecre ase in d compared to 
d/b = \/l37/4 as used in equation 
similar to taking C ■ 
for our purposes. 

We comment that a more general torque formula, valid 
for higher masses, would therefore have C = C(b/h,q) and 
d = d(q). A detailed analysis is beyond the scope of this 
paper, but we point out that the results of|Masset et al.| 
(2006) and Paardekooper & Papaloizou (2009b I could be 
used to derive approximate functional forms for the mass 
dependence of C and d. 



7 DISCUSSION 

In equation |47| we have presented a formula for the un- 
saturated torque on low-mass planets. This torque consists 
of the Lindblad torque, together with the barotropic and 
entropy-related horseshoe drag. The barotropic part of the 
horseshoe drag is due to material conserving its vortensity, 
and its expression is identical to the barotropic case. The 
entropy-related part of the horseshoe drag is exerted by den- 
sity structures produced by material conserving its entropy, 
plus an additional component linked to the production of 
specific vorticity at the outgoing separatrices. Masset & Ca- 
soli ( 2009 1 recently studied the case with (3 — in detail, 



arguing that there is no contribution to the torque associ- 
ated with a density response resulting from entropy advec- 
tion. If we neglect this term in our simple model, and take 
v — 3/2, our expression for the density perturbation agrees 
with that in Masset & Casoli (20091. We comment that the 



contribution of entropy advection to the total torque is small 
compared to that due to vortensity generation (20 — 30%, de- 
pending on the softening length), but keeping it gives better 
agreement with numerical simulations. 

Since we consider the unsaturated torque, equation |47| 
should apply only in regions of the discs where thermal and 
viscous diffusion keep the corotation torque unsaturated. We 
will study saturation effects in a forthcoming work. Here, we 
just comment that it has been shown ([Paardeko oper fc Pa-| 
paloizou|2008 Kley & Crida 2008 ) that when including vis- 
cosity as well as thermal diffusion (possibly through radia- 
tive effects) a sizable fraction of the unsaturated corotation 
torque can be sustained. 

We have worked in the 2D approximation throughout 
this paper. Although we have used a reasonable value for 
the gravitational softening parameter to mimic 3D averag- 
ing, a fully 3D model of the horseshoe region is required to 
capture possible effects due to vertical motions. The strong 
dependence of the torque on softening suggests that non- 
isothermal effects in 3D may be very strong, but it is im- 
portant to keep in mind that the torque depends on the 
detailed flow structure around the planet, which may well 
be different in 3D. 

We have neglected effects of any magnetic fields. It re- 
mains to be seen, for example, whether a fully turbulent disc 



in Kley & Crida (20081 



( Nelson & Papaloizou 2004 1 allows for horseshoe turns to oc- 
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cur. Self-gravity was also ignored. It was shown in |Pierens] 
& Hure ( 2005 1 that self-gravity tends to make the wave 



torque slightly stronger due to a shift in the Lindblad reso- 
nances. This was confirmed numerically by[Baruteau fe Mas-| 



set 



(2008b), who also showed that the impact of self-gravity 



on the corotation torque is small. 

When calculating the torque on the planet, we have in- 
cluded all disc material. Tests have shown that it makes very 
little difference for these low-mass planets when a fraction 
of the Hill sphere is excluded. The situation is different for 
high-mass planets, for which a circumplanetary disc may 
appear. In those cases it is an issue which material should 
exert a torque on the planet 1 Crida et al.||200"9 l. However, 
even the highest mass planets we consider do not show any 
evidence for a circumplanetary disc around the planet (see 



the bottom panels of Fig. 131, so this issue is of no concern 
here. 

We have kept the planet on a fixed circular orbit. There- 
fore, we have neglected any distortion of the streamline 
topology due to the radial movement of the planet. This 
can have some impact on the corotation torque, especially 
for massive discs in which the planet migrates fast enough, 
under influence of the Type I torque discussed in this paper, 
so that r p ~ x s Q p . Then one may expect to see migration 
behaviour similar to Type III. This has not been considered 
so far. 

It is important to note that Type I migration will al- 
ways be fast, unless the background disc is close to the zero- 
torque lines in Fig. [7| The torque predicted by equation [47] 
can indeed be much larger in magnitude than the linear, 



isothermal Type I torque (see equation 20 1 . Type I migra- 



tion can be directed inward or outward, depending on the 
background entropy gradient. Outward migration is always 
limited, however, since inevitably the planet will enter a re- 
gion of the disc where the opacity is low enough to make 
cooling efficient, pushing the planet back into the (locally) 
isothermal regime of inward migration. 

One scenario that permits slow migration only is the fol- 
lowing. If the thermodynamic state of the inner disc is such 
that it permits outward Type I migration, there exists an 



equilibrium radius r e where the torque is zero ( Paardekooper 
fc~Me llcma 2008). A low -mass planet will then tend to mi- 
grate towards r e , either from the outer disc or from the inner 
disc. This radius r c will move inward when the disc is losing 
mass, either by accretion onto the star or by evaporation, 
taking the planet along. This way, low-mass planets can mi- 
grate slowly (on a time scale comparable to the disc life 
time) towards the central star. 



8 CONCLUSIONS 



We have presented a simple relation (equation 47 1 that gov- 
erns the migration speed and direction for low-mass planets. 
Since we have considered unsaturated torques only, this law 
should apply in regions of the disc where thermal and vis- 
cous diffusion act to keep the corotation torque unsaturated. 
The total torque is found to strongly depend on the presence 
of a radial entropy gradient in the disc, with the possibil- 
ity of outward migration in the case of outward decreasing 
entropy. 
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